/*
 *  MWDots.cpp
 *  T3nsors2
 *
 *  Created by Michael Barriault on 10-06-30.
 *  Copyright 2010 University of Guelph. All rights reserved.
 *
 */

#include <cmath>
#include "MW.h"
using namespace t3;

scalar* MW::PhiDot(scalar* Phi, scalar* Psi, vector* E, vector* B) {
	scalar* PhiDot = scalar::init(O, "PhiDot");
	scalar* dE = div(E);
	scalar* Phi4 = lapfour(Phi);
	FOR(o,O->size()) PhiDot->at(o) = -K*Phi->at(o) + dE->at(o) + 4*M_PI*rho->at(o) - epsilon*pow(O->d(0),3.)*Phi4->at(o);
	dE->release(); Phi4->release();
	return PhiDot;
}

scalar* MW::PsiDot(scalar* Phi, scalar* Psi, vector* E, vector* B) {
	scalar* PsiDot = scalar::init(O, "PsiDot");
	scalar* dB = div(B);
	scalar* Psi4 = lapfour(Psi);
	FOR(o,O->size()) PsiDot->at(o) = -K*Psi->at(o) + dB->at(o) - epsilon*pow(O->d(0),3.)*Psi4->at(o);
	dB->release(); Psi4->release();
	return PsiDot;
}

vector* MW::EDot(scalar* Phi, scalar* Psi, vector* E, vector* B) {
	vector* EDot = vector::init(O, 3, "EDot");
	vector* dB = curl(B);
	vector* dPhi = diff(Phi);
	vector* E4 = lapfour(E);
	FOR(a,3) FOR(o,O->size()) EDot->c(a)->at(o) = dB->c(a)->at(o) + dPhi->c(a)->at(o) - 4*M_PI*J->c(a)->at(o) - epsilon*pow(O->d(0),3.)*E4->c(a)->at(o);
	dB->release(); dPhi->release(); E4->release();
	return EDot;
}

vector* MW::BDot(scalar* Phi, scalar* Psi, vector* E, vector* B) {
	vector* BDot = vector::init(O, 3, "BDot");
	vector* dE = curl(E);
	vector* dPsi = diff(Psi);
	vector* B4 = lapfour(B);
	FOR(a,3) FOR(o,O->size()) BDot->c(a)->at(o) = -dE->c(a)->at(o) + dPsi->c(a)->at(o) - epsilon*pow(O->d(0),3.)*B4->c(a)->at(o);
	dE->release(); dPsi->release(); B4->release();
	return BDot;
}

